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Abstract 

We extend the notion of numerical stability of finite difference approximations to 
include hyperbolic systems that are first order in time and second order in space, 
such as those that appear in Numerical Relativity and, more generally, in Hamil¬ 
tonian formulations of field theories. By analyzing the symbol of the second order 
system, we obtain necessary and sufficient conditions for stability in a discrete norm 
containing one-sided difference operators. We prove stability for certain toy models 
and the linearized Nagy-Ortiz-Reula formulation of Einstein’s equations. 

We also find that, unlike in the fully first order case, standard discretizations of 
some well-posed problems lead to unstable schemes and that the Courant limits are 
not always simply related to the characteristic speeds of the continuum problem. 
Finally, we propose methods for testing stability for second order in space hyperbolic 
systems. 
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1 Introduction 


The Einstein equations consist of a set of ten coupled non linear second order 
partial differential equations. In order to perform numerical time evolutions 
the fully second order system is usually written as a first order in time sys¬ 
tem. Such systems can be evolved directly [1,2], or a further reduction from 
second to first spatial order can be performed (see, for example, [3,4,5,6]). 
Whereas the theory of Cauchy problems for fully first order systems of partial 
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differential equations is understood, in terms of well-posedness at the con¬ 
tinuum and the stability of finite difference approximations, the theory of 
second order in space hyperbolic systems is less well developed. The recent 
improvement in the understanding of second order in space formulations of 
Einstein’s equations at the continuum [7,8,9,10,11], has not been matched by 
developments concerning finite difference approximations of such systems (see, 
however, [12,13]). Given that these systems have fewer variables, fewer con¬ 
straints, and typically smaller errors (see [12] and Appendix B), it is desirable 
to better appreciate their properties. Note that first order in time hyperbolic 
systems, which are not necessarily Erst order in space, also arise naturally in 
Hamiltonian formulations of held theories. 

The standard notion of stability for fully first order systems based on the 
discrete L 2 norm is unsuitable for analyzing second order in space hyperbolic 
systems. This can be understood by analogy with the continuum result for the 
one dimensional wave equation written in first order in time and second order 
in space form: d t </>(t, x) = U(t, x ), d t Tl(t, x) = <9]i</>(f, x). Consider the family of 
solutions 4>(x,t) = sin(o;x) cos(o;t), tt( x,t) = — a;sin(a;x) sin(uA) generated by 
the initial data 4>o(x) = sin(a;x), 7io(x) = 0. By varying u) in the initial data, 
the L 2 norm of the solution at a fixed time t, / 0 27r ( \4>\ 2 + |n| J )<ia;, can be made 
arbitrarily large with respect to the initial data (whose norm is independent 
of to), thus contradicting well-posedness of the Cauchy problem in L 2 [14,15]. 
The introduction of the new variable, X = d x (j) : allows the construction of 
a first order system, the Cauchy problem of which is well-posed in L 2 . The 
original second order problem can then be shown to be well-posed in a norm 
containing derivatives, namely Jo 27r (|0| J + |Id| 2 + \d x (j)\ 2 )dx, which corresponds 
to the L 2 norm of the first order reduction. 

In this work we consider linear constant coefficient Cauchy problems. We use 
the method of lines to separate the time integration from the spatial discretiza¬ 
tion. We show that by reducing the discrete system to first order in Fourier 
space, it is possible to determine stability in physical space with respect to 
a discrete norm containing one-sided difference operators. This is done by 
extending the notion of a symmetrizer to the discrete case. We apply these 
techniques to problems, starting with the wave equation written as a first 
order in time, second order in space system. We consider both second and 
fourth order accurate discretizations. A similar but more complicated analysis 
is done for the Knapp-Walker-Baumgarte (KWB) [16] and Z1 [17] formula¬ 
tions of electromagnetism, and the Nagy-Ortiz-Reula (NOR) [8] formulation 
of Einstein’s equations. We also point out stability issues related to the ADM 
[18] and Z4 [19] formulations. 

In Sec. 2, we summarize some relevant material from the literature. In Sec. 3 
we introduce the concept of a discrete symmetrizer. We also illustrate the 
reduction procedure to first order in Fourier space, which can be used for 
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obtaining energy estimates at the continuum. We introduce the analogous 
idea for the discrete case, and discuss convergence. In Sec. 4 we apply these 
techniques to the systems mentioned above. We propose methods in Sec. 5 
for testing stability experimentally both for linear and non linear systems. We 
summarize the main results of this paper in Sec. 6. In Appendix A, we describe 
the different time integration methods that we consider, and in Appendix B 
we compare numerical properties of the wave equation written as a first order 
system with those of the wave equation written as a first order in time, second 
order in space system. In Appendix C we highlight differences in the constraint 
propagation properties between first and second order systems. 


2 Background 


Well-posedness, the (local in time) existence of a unique solution which de¬ 
pends continuously on the problem’s data, is a fundamental requirement for 
the successful generation of numerical solutions approximating the solution of 
a continuum problem. In this section we review the notion of well-posedness for 
linear constant coefficient Cauchy problems, as well as the concept of stability 
for finite difference approximations. We conclude the section by providing a 
simple sufficient condition for stability of first order fully discrete problems 
based on the properties of the symbol of the semi-discrete system, which will 
be extended to discretizations of second order in space problems in the next 
section. 


2.1 Constant coefficient Cauchy problems 


In this work we will be dealing with initial value (or Cauchy) problems of the 
form 


^, X )= p (cy {ttXh (1) 

MO ,x) = f{x), (2) 

in d spatial dimensions, where x G R d , u = (u^\u^ 2 \ .. .u^) T and P is a 
linear, constant coefficient, differential operator of order p. We consider only 
the cases p — 1 and p = 2. Furthermore, we assume that the eigenvalues of 
the symbol of the differential operator, P{iuj ), which is obtained by replac¬ 
ing d/dxj in P{d/dx) with iujj, for j = 1,2,... ,d, have real part uniformly 
bounded from below and above. We are thus excluding parabolic systems, but 
we are allowing for systems like the wave equation written as a first order in 
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time, second order in space system. For simplicity we focus on solutions that 
are 27r-periodic in all spatial coordinate directions. Thus the initial data, f(x), 
is chosen so that it satisfies this property. 

We consider the p = 1 case, leaving the p — 2 case for the next section. 
Following Definition 4.1.1 in [20] we say that problem (l)-(2) is well-posed with 
respect to a norm || ■ || if for every smooth periodic / there is a unique smooth 
spatially periodic solution and there are constants a and K , independent of 
/, such that for t > 0 

\\u(t,-)\\<Ke at \\f\\. (3) 


Exponential growth must be allowed if one wants to treat problems with 
lower order terms. For first order hyperbolic systems the L 2 norm ||u>|| 2 = 
/ 0 27r ... / 0 27r \w{x)\ 2 dx\... dxd is usually used in (3). We will see later that the 
second order systems we study in this work require the use of a different norm. 

Taking /(x) = ( 2n )~ d / 2 f(uj) the formal solution of (l)-(2) is u(t, x) = 

(2n)~ d ^ 2 J2uJ el< ' w ’ x ' >e ^ u> ^ t f( UJ )- K can be shown (Theorem 4.5.1 in [20]) that 
well-posedness in the L 2 norm is equivalent to there being constants K , a 
such that, for all uj and for t > 0, 

| e PMt| < Ke at ; (4) 


where \A\ — sup| u | =1 \Au\ is the matrix (operator) norm of a matrix A. 

Well-posedness of the Cauchy problem in the L 2 norm is also equivalent (The¬ 
orem 4.5.8 in [20]) to the existence of constants a, K > 0 and of Hermitian 
matrices H (cn) satisfying 1 , for every u, 


K~ l I < H(u) < KI , (5) 

H{uj)P{iuj) + P*(iu))H(u>) < 2 aH(u ), 

where P* represents the Hermitian conjugate of P. The last inequality gives 
an energy estimate for each Fourier mode and the estimate in physical space, 
Eq. (3), follows from Parseval’s relation, \\u(t, -)|| 2 = |w(f> a; )| 2 - Since the 

existence of H(u>) is not affected by the addition of a constant matrix to 
P(iu>) (Lemma 2.3.5 in [21]), undifferentiated terms on the right hand side of 
the equations can be ignored in the analysis of well-posedness. If (5) is satisfied 
with HP + P*H = 0 then H is called a symmetrizer. 

1 Two Hermitian matrices, A and B, satisfy A < B if and only if x*Ax < x*Bx 
for every x. If a matrix H(uj) satisfies K~ l I < H(oj) < KI for every ui, we say that 
H(lo) is equivalent to the identity matrix. 
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For p = 1, system (1) is said to be strongly hyperbolic if the correspond¬ 
ing Cauchy problem is well-posed in the L 2 norm (i.e. if H{ui) exists) 2 . If 
H{u>) = /, the system is said to be symmetric hyperbolic. If H{u) — H is 
independent of uj, then we say that the system is symmetrizable hyperbolic 3 . 

In this case the change of variables u = H 1 ^ 2 u brings the system into symmet¬ 
ric hyperbolic form. Finally, well-posedness is not affected by the presence of 
forcing (inhomogeneous) terms (Theorem 4.7.2 in [20]). For cases where such 
terms are present, the estimate requires modification. 

Note that, in the absence of lower order terms, whereas symmetrizable hy- 
perbolicity guarantees the existence of a conserved energy in physical space, 
(u, Hu), a strongly hyperbolic system satisfies the estimate \\u(t ,-)|| < K\\u(0 ,-)|| 
with a constant K > 1. Furthermore, in the variable coefficient case, well- 
posedness results require smoothness of the synnnetrizer H(x,t,cu) in all ar¬ 
guments [21]. 


2.2 Numerical stability 

2.2.1 Notation 

Our notation and conventions follow closely those of [20]. We introduce a spa¬ 
tial grid Xj = .. .,xfj) = {jih l ,j 2 h 2 , ■ ■ ■ ,jdh d ), where h r = 2vr/iV r 

and j r = 0,1,..., N r — 1, and the vector-valued grid function Vj(t) approxi¬ 
mating u(t, Xj). Periodicity requires that Vj = v mo d(j,N)- The partial derivatives 
in (1) are approximated using either the standard second order accurate dis¬ 
cretization 


a a 


0 i ' 


| D 0i D 0j if i a .1 
\ D+iD^i if i = j 


or the standard fourth order accurate discretization 


( 6 ) 


8, - D< 4 > = (1 - y D *> D ~) ' (7) 

2 For dtu = A l diU the symbol is P = iojiA l . The system is said to be weakly hyper¬ 
bolic if the eigenvalues of P{iuj) are imaginary. Strong hyperbolicity is equivalent 
to P{iuj ) being uniformly diagonalizable with imaginary eigenvalues. We define the 
characteristic speeds in the direction uji to be the eigenvalues of P(iu>) divided by 
iuj. 

3 Symmetrizable hyperbolic systems are often also called symmetric hyperbolic. 


5 



didj 


' DfDf if i+j 

< 

D+iD-i (l - £D +l D^) if i = j 


J 


where D + Vj = (vj+ 1 — Vj)/h, D_Vj = (vj — Vj-i)/h, D 0 Vj = (vj+ 1 — Vj-i)/2h , 
and D + D_Vj = (u J+1 — 2u,- + Vj_i)/h 2 . The discretization of <9 2 as in (6) or 
(7) gives the desired order of local accuracy without requiring a larger sten¬ 
cil. We then integrate the resulting system of rri IXLi N r ordinary differential 
equations 


j t v j {t) = Pv j {t), (8) 

Vj(0) = fj, ( 9 ) 

where fj = f{xf), with three different time integrators. These are iterative 
Crank Nicholson (ICN) and third and fourth order Runge-Kutta (3RK and 
4RK) methods, which are widely used by numerical relativists (see Appendix 
A for definitions). Using the fact that the operator P is linear and time in¬ 
dependent we can write the fully discrete system in polynomial form (see for 
example [20]) 


vr = Qv] = V(kP)v h ( 10 ) 

«? = /;, ( 11 ) 

where k = Xh is the time step, A is called the Courant factor, and u” represents 
the grid-function at time t n = nk. This is an explicit, one step, scheme. For 
ICN we have V(x) = 1 + 2§f, whereas for p -th order Runge-Kutta we 
have V{x) = £r=o fr- 


2.2.2 Definition of stability 

We recall the definition of numerical stability and discuss some necessary and 
sufficient conditions. The solution of the finite difference scheme (10)—(11) is 
v n = Q n f. We introduce the scalar product ( u,v)h = l>2j{uj,Vj)h d , where 
h d = nil hi, j = (ji, j 2 , ■ ■ ■ ,jd) is a multi-index and {uj,vf) = upvf'f 
This allows us to define a norm ||u||/i = (v,v)]/ 2 . The approximation (10)—(11) 
is said to be stable with respect to this norm if there exist constants a, K, 
such that for all h, k, 0 < h < h 0 , 0 < k < k 0 , the estimate 

\\v n \\ h <Ke at "\\f\\ h (12) 

holds for all n such that t n = nk and all initial grid-functions /. This concept 
of stability is the discrete analogue of (3). It guarantees that the solutions are 
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bounded as h —> 0. However, the schemes we consider are at most conditionally 
stable. By this we mean that there exists a Ao such that the above inequality 
holds if and only if the additional condition A = k/h < A 0 is satisfied. 

Theorem 5.1.2 in [20] guarantees that if the scheme (10)-(11) is stable, then 
the modified scheme 


< +1 = {Q + kR)v™ , (13) 

«? = /; ( 14 ) 

is also stable provided that R is bounded. This will be the case when R rep¬ 
resents constant terms (lower order terms) in the continuum problem. Hence 
for a first order in space system lower order terms can be ignored without 
affecting stability. 


2.2.3 Convergence 

Following Theorem 5.1.3 in [20], consistency and stability imply convergence. 
Assume that the continuum solution u of (l)-(2) is smooth and that the 
scheme (10)-(11) is stable. Further assume that the scheme and the initial 
data are consistent. Then, on any finite interval [0, T ], the error satisfies 

\\v n -u(;t n )\\ h <0(h P '+k P *) (15) 


i.e. the solutions of the finite difference scheme converge as h —> 0 to the 
solution of the differential equation 4 . 


2.2.4 Fourier analysis of stability 


For approximations with constant coefficients, Fourier analysis can be used 
to obtain necessary and sufficient conditions for stability which can be more 
easily verified than the above definition. We assume that N, the number of 
grid-points in each direction, is even (the odd case is discussed in Sec. 2.2.5). 
If we represent vf by 


1 

(2-7t)^ 


5>^>u>), 


(16) 


4 Note that the big O in inequality (15) contains higher derivatives of the exact 
solution. Smoothness of the solution of the continuum problem is not required for 
convergence. For instance, a weaker condition for fourth order convergence (p± = 
P 2 = 4) is that the solution be C 5 . 
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where uj = ( 001 , 002 , ■ ■ ■ ,u>d), oo r = —N/2 + 1,..., N/2, and substitute it into 
the difference scheme (10)—(11), we obtain 


«"+V) = Q(0*», (17) 

«Vo=/M, (is) 

where = oo r h = —n + 2tt/N, —it + Ati/N, ..., +7r and r = 1, 2,..., d. The 
m x m matrix 0(0 is called the amplification matrix of the scheme and is 
a real polynomial in P, the symbol of the Fourier transformed semi-discrete 
problem, 

0(0 = P(kP(£)). (19) 


The matrix P(£) will play an important role in the next section. It can be 
readily computed from P in Eq. (8) with the replacements 


1 

Doi T sin 0, 
h 

(20) 

n n ^ • 2 0 

D +i D_i — > — —r sin —. 

+1 1 h 2 2 

(21) 

Using the discrete Parseval’s relation 

IHI 2 , = EI^MI 2 

UJ 

(22) 


and the fact that the solution of (17)-(18) is v n (u) = Q n (£)f(ou) one can show 
(Theorem 5.2.1 of [20]) that a necessary and sufficient condition for stability 
with respect to the || • [[^ norm is given by 

\Q n (0\<Ke atn (23) 

for all h = 2n/N < ho, k < fc 0 , n with t n = nk, and ui r = — N/2 + 1,..., N/2, 
r = 1,2,... ,d. 

A much easier condition to verify is the von Neumann condition, which is only 
a necessary condition for stability. It corresponds to the requirement that the 
eigenvalues z v (£) of 0(0 satisfy 

MOI<e afe (24) 


for all h < ho and |£ r | < n. However, when the amplihcation matrix can 
be uniformly diagonalized (i.e. there exists a non-singular matrix T(£) that 


diagonalizes Q(£) and satisfies |T(£)| |T _1 (£)| < C with C independent of 0 
then the von Neumann condition is also sufficient for stability. In particular, 
if Q is normal then it can be unitarily (and therefore uniformly) diagonalized, 
|T(0| = |r _1 (0l = 1- Since for the time integrators that we consider Q is a 
polynomial in P, Q will be normal if P is normal (as would be the case if P 
were Hermitian or anti-Hermitian). Note that if the von Neumann condition 
is violated then the scheme is not stable in any sense. 

It is possible for a discretization to be (conditionally) stable without Q being 
normal (and hence unitarily diagonalizable). This turns out to be the case for 
most systems considered in this work. In such cases we fold it convenient to 
introduce the norm \u\jj = {u 1 Hu) 1 ^ 2 * * and proceed as follows. Let us assume 
that //”(£) are Hermitian matrices such that 


K~ X I < H(£) < KI, (25) 

\Q\h < e* 


where K is a positive constant. Notice that 5 \u\fy = \H 1 ^ 2 u\ < K 1/2 \u\ and 
K~ l \A\ < \A\fj = \H 1/,2 AH~ 1 / 2 \ < A'|A|. As a consequence the von Neumann 
condition is satisfied, a(Q) = a^H^-^QH" 1 ^ 2 ) < \H 1 / 2 QH~ 1 / 2 \ — \Q\^ < e ak , 
where a(Q ) denotes the spectral radius of Q. Stability follows from 

\Q n \ < K\Q n \ k < K\Q\ n k < Ke at (26) 


According to the Kreiss Matrix Theorem (Sec. 4.9 of [14]), for a family T of 
m x m matrices A the following two statements are equivalent: 

(1) There exists a constant C such that for all A e T and all positive integers 
n 


\A n \ < C 

(2) There is a constant K > 0 and, for each A G A, a positive definite 
Hermitian matrix H with the properties 

A" -1 / < H < KI, A*HA < H 

This implies that condition (26) is also necessary for stability. 


5 For a positive definite Hermitian matrix H, H a (for a not necessarily an integer) 
is defined as S*D a S where H = S*DS and D is the diagonal matrix of positive real 
eigenvalues. 
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2.2.5 Number of grid points 


In this review we have assumed that the number of grid points in each direction 
is even. This means that no matter how small the number of grid points is, as 
long as it is even, the highest frequency £ r = n is present. To allow for an odd 
number of grid points one must change the summation range in Eq. (16) to 
ay = — (N — l)/2,... (N — l)/2, in which case, |£ r | never equals ir, although 
it does approach this value as h —> 0. 


2.3 A sufficient condition for stability 


We can now give a simpler sufficient condition for numerical stability. This 
condition applies to systems which admit a conserved energy in Fourier space 
and will enable us in Sec. 3.2 to obtain another condition suitable for the 
applications. We consider only time integrators such that 

Q = VikP) (27) 


The eigenvalues q v of Q are related to the eigenvalues p u of P by q u = V{kp v ). 
This can be seen by using Shur’s lemma. Provided that the eigenvalues p u are 
imaginary, the inequality \q v \ < 1 is equivalent to kp v < ao, where «o = 2 for 
ICN, for 4RK, for 3RK. Hence, 


a(kP) < «o 


(28) 


is equivalent to cr(Q) < 1. Condition (28) is called local stability on the imag¬ 
inary axis in [22], Suppose that the time step is such that a(kP) < a 0 . If we 
can End Hermitian matrices such that 


K~ l I < H{£) < KI , (29) 

h(0P(0 + P(0*h(0 = o, (30) 

we say that H(£) is a discrete symmetrizer of P(£). The matrices H X ^ 2 PH V ^ 2 
are anti-Hermitian, hence they can be diagonalized by unitary matrices S(£). 
This implies that the matrices diagonalize Q(£). The inequality 


\Q\h = \H 1/2 QH~ 1/2 \ = | S- l H l/2 Qpl- ll2 S\ = <j(Q) < 1 (31) 

guarantees stability. In fact, the amplification matrix can be uniformly diag¬ 
onalized by T(0 = P" 1/2 (£)S(£). 
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In applications one would construct a norm (i.e., matrices H(£) satisfying 
(29)) which is conserved by the Fourier transformed semi-discrete evolution 
equations, 


J t \v\\ = («, (HP+P’H) 6)=0. 


(32) 


This implies that condition (30) holds and H{fi) is a discrete symmetrizer. 

To construct H one can proceed as follows. Assume the existence of a matrix 
T such that T~ 1 PT = A is diagonal with imaginary elements. Then the quan¬ 
tity v*Hv, where H = T” 1 * DT^ 1 and D is a positive definite matrix which 
commutes with A, is conserved by the system dtv = Pv. Defining the char¬ 
acteristic variables of P to be w = T~ l v (these are individually conserved: 
d t \wi\ 2 = 0), we see that to construct a conserved quantity one can take w*Dw. 
(For D = I this corresponds to adding the squared absolute values of the char¬ 
acteristic variables.) For H to be a symmetrizer it remains to be established 
that K~ l \v\ 2 < v*Hv < K\v\ 2 . 


3 Stability of first order in time, second order in space systems 


What we have done so far applies to fully first order systems. We have shown 
that if inequalities (28) and (29) and Eq. (32) hold, then the fully discrete 
scheme is stable and satisfies the estimate (12) with a = 0. In this section we 
show how this can be extended to second order in space systems. We first look 
at the continuum problem and then investigate its standard discretization. 


3.1 Well-posedness of first order in time and second order in space hyperbolic 
systems 


It is possible for the Cauchy problem of a first order in time and second order 
in space system of equations to be ill-posed in the L 2 norm, but well-posed 
in a norm which contains additional derivatives (see the introduction). The 
system is still useful; for example, a suitable finite difference approximation 
of the equations can be convergent in the discrete L 2 norm. We analyze the 
well-posedness of the Cauchy problem for such systems by using the analytical 
tool of a reduction to first order. This will be done in Fourier space, so that 
the number of additional variables being introduced is minimized [23]. 

Consider system (1) with p = 2 and suppose that it can be written in the form 
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d t u = Pu , 


u = 


C 

G l di + J 


(33) 


P 



1 A% + B 

y D lj didj + E l di + F 



where the evolved variables have been split into two types. The column vector 
u represents those that are differentiated twice (in space) and v represents 
those that are not. In P a sum over repeated indices is assumed. Not all second 
order in space systems can be written in this form (for example, u t = u xx ). 
This form is general enough to include all the first order in time, second order 
in space systems that we have considered that can be reduced to first order in 
space. Fourier transforming this system, we obtain 


d t ii 


P 


Pu. 



' iujA n + B 

i —u 2 D nn + iujE n + F 


C 


iu)G n + J 


(34) 


where M n = and a= \uj\rii and uj = |u;|. We define the second order 

principal symbol to be 


P' 


1 icuA n 
y -u 2 D nn 



(35) 


We now state the main result of this subsection. If there exists H{uj) = H*(uj ) 
such that the energy u*Hii is conserved by the principal system d t u = P'ii 
and H satisfies 


K- X F <H< KI m F 


V° 



(36) 


where K is a positive scalar constant, then the solution of (33) satisfies the 
estimate 


||«(V)II < Ke at \\u(0, -)||, 


(37) 
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||w|| 2 = / |w| J + \dju\ 2 + \v\ 2 d d x , 

J i =1 

and the problem is well-posed in this norm 6 . 

The proof proceeds via a pseudo-differential reduction to first order [8]. This 
involves the introduction of a new variable w = icon. By taking a time deriva¬ 
tive of this definition, we obtain the enlarged system in which the second 
derivative of u has been replaced with a first derivative of w. We reduce the 
order of the system as much as possible so that any occurrence of icuu is 
replaced with w. This particular first order reduction is 


d t u R 


PrUr , 


Ur = 




W 


Pr — 


' B A n C N 

0 icuA n + B iuC 
F iujD nn + E n iu)G n + J ^ 


(38) 


This system is equivalent to the second order system (34) only when the 
auxiliary constraints 

C(t,u) = w(t,ui) — iuu(t,u) = 0 (39) 


are satisfied. It can be shown that d t C = BC so if these constraints are satisfied 
initially, then they are satisfied for all time. They are said to be propagated by 
the first order evolution equations. 

If this system admits a matrix Hr satisfying (5) then the solutions satisfy the 
estimates 


\iiR{t,uj)\ < Ke at \u R (0,uu)\ , (40) 

where \ur \ 2 = |u| 2 + |w)j J + |u| 2 , for arbitrary initial data and u. Specifically, 
the estimate holds for solutions which satisfy the auxiliary constraints and 

6 Note that we made no assumptions regarding the smoothness of the matrix H(uj). 
In view of generalizations of this work to the variable coefficient case it may be 
desirable to demand that T -1 * HT ^ 1 , where T is defined in Eq. (44), be smooth in 
all variables. 
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therefore correspond to solutions of the second order system. The uniform 
estimate in uj of 


\u \ 2 + uj 2 \u \ 2 + | u | 2 = |«| 2 + ^2 K < W &| 2 + l ^| 2 ( 41 ) 

i —1 

implies, by Parseval’s relation, the estimate in real space 


JNV)|| <^e°*||'u(0,-)||, (42) 

m| j + \diu\ 2 + \v\ 2 d d x. 

2—1 

So the existence of Hr for a first order pseudo-differential reduction implies the 
well-posedness of the second order system with respect to a norm containing 
derivatives. 

We have still to show that we can find an Hr for (38). Whether or not this 
is the case is independent of the lower order terms Pr contains. A calculation 
similar to Lemma 2.3.5 in [21] shows that if P{eo) admits an Hr , then so will 
P{ui) + B(uj), where B{uo) is any matrix which satisfies \B\ + \B*\ < C for C 
independent of uj. In other words, the terms that are not multiplied by iu can 
be dropped from (38), giving the principal symbol of the first order reduction 



Pr 


0 0 0 
0 iujA n icuC 
0 iuD nn iujG r 


(43) 


without affecting the well-posedness. The principal symbols of the second order 
system, Eq. (35), and the first order pseudo-differential reduction, Eq. (43), 
are related by 


P'r = 


^0 

0 \ 


* iuj 0 \ 



- 4 h 

T = 

• 

(44) 

v° 

TP'T~ l J 


1° V 



(Note that T 1 does not exist for uj = 0. However, in this case, Pr = 0, and 
admits the identity as a symmetrizer.) By assumption, there exists H(uj) = 
H*{uo) such that u*Hii is conserved by the principal system d t u = P'u and 
satisfies (36). This H satisfies HP' + P'*H = 0, and it is straightforward to 
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show that 


0 \ 

(45) 

T-^HT- 1 j 

satisfies H R = H R and H R P' R + 77# = 0. Further, by noting that T*T = 7^, 

using (36) one can show^ that H R satisfies 77 7 ^ Hr < KI. Hence we have 
found a symmetrizer of and the result has been proved ' . 

To construct H one can use the characteristic variables of P', as described 
at the end of Sec. 2.3. We would like to point out that this analysis did not 
require that the auxiliary constraint propagation problem be well-posed. These 
constraints are merely a tool for the analysis of the system. We only need to 
establish uniqueness of the solution with zero initial data for the auxiliary 
constraints. In the linear constant coefficient case this result is trivial. When 
evolving the second order system, these constraints are identically zero at all 
times. An alternative to the pseudo-differential reduction method is to perform 
a fully differential reduction by introducing a new variable in physical space 
for each derivative (see, for example, [7,11]). 


Hr = 


3.2 Stability of discretizations of first order in time and second order in space 
systems 


We now show how the continuum analysis of the previous subsection can be 
extended to the fully discrete case. The semi-discrete finite difference approx¬ 
imation of (33) can be written as 


d 

dt 


v = Pv , 



f A* D*P + B 
y I) ;i I))f + + F 


G i D ( f > + JJ 


(46) 


where D \ 1 ^ is a discretization of the first derivative in the i direction and 
D)fi is a discretization of the second derivative in the i and j directions. For 
example, the standard second order accurate discretization would have 

' It can also be shown that Pr is diagonalizable with the same eigenvalues as P', 
plus as many zeroes as there are components of u. 
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dP = a 


Oii 


d (2) = 

l] 


DoiD 0 j i j 

D +1 D_ t i = j 

The principal symbol of the semi-discrete system is 


(47) 


P' = 


( At)\ l) 

I)'A)A 


c 


G l D 


(i) 


(48) 


where 


A W = ^sin&, Dif = 


h 2 
4 

7^ 


sin & sin 6 i 7 ^ j 


• 2 Si 
Sill — 


(49) 


t = J 


for the standard second order discretization. The pseudo-discrete first order 
reduction is obtained by defining 


w = Alii , A = \D 


+i | 


Z =1 


The reduced system is 


(50) 


-rVR = PrVr, 

dt 


v R = 




w 

\VJ 


Pr = 


-1 Ai n ( 1 ) 


B (ity-'AD 
0 APP + B 


(51) 


\^F (inypD'Wff + E'D\ L) ) g i d\ l> + j) 


c 

ific 

E l Dp ] ) G l D ( P 


(52) 


We can show that the discrete auxiliary constraint is preserved by the time 
integrator. Define c = (—iQI I 0), so that the constraint is cvr = w — iVtu = 0. 
Since cPrVr = Bcvr, we have that cvr = 0 implies cPrVr = 0 and hence 
cP r Vr = 0 and cP(kP)vR = 0. Now consider evolving the reduced system 
with a polynomial time integrator; i.e. u ^ +1 = V{kpR)v n R . If the auxiliary 
constraints are satisfied on one time step, then they are satisfied on the next 
as well, since cv n R = 0 implies cv 1 ^ 1 = cV(kP)‘v n R = 0 . Hence there is a one-to- 
one correspondence between solutions of the second order fully discrete system 
and those of the constraint-satisfying reduced system. Note that we have used 
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the fact that the time integrator is a polynomial in Pr, as is the case for 
systems with constant coefficients. This result can be extended to the variable 
coefficient case, where one would have to perform the reduction to first order 
in physical space by introducing the gridfunctions = D +i u. 

Making use of Theorem 5.1.2 of [20], the terms which correspond to the contin¬ 
uum lower order terms can be dropped from Pr without affecting the stability 
of the fully discrete system, provided that (zh 2 )~ 1 P>j 1 \ kt)^ } and are 

bounded. This guarantees that the assumptions of the theorem are satisfied. 
This is true for the second and fourth order accurate standard discretizations. 

The result for stability of the fully discrete problem is analogous to that for 
well-posedness at the continuum. If there exists P(£) = P*(£) such that the 
energy v*Hv is conserved by the semi-discrete principal system d f v = P'v 
and H satisfies 


K - 1 ! n < H < KI n , 


In = 


/ n 2 

V° 



(53) 


where K is a positive scalar constant, then it is possible to construct a discrete 
symmetrizer for the first order reduction with no lower order terms. So if, in 
addition, the principal symbol P' satisfies a(kP') < ck 0 ; the fully discrete 
system (including lower order terms) is stable with respect to the norm 


HI h,D + = Mil + Ml + \\ D +i u \\h, 

i=l 

i.e. the solution satisfies the estimate 


(54) 


|t> n |k, Di < KC 


OLtn 


V 


\h,D^ 


(55) 


Again, H can be constructed from the characteristic variables of P', as de¬ 
scribed at the end of Sec. 2.3. Note that the matrix Pr is not defined for 
fl = 0. However, this does not cause any difficulties in the linear constant co¬ 
efficient case. One can write the space of solutions as a direct sum consisting 
of constant functions plus a space of solutions with nontrivial O, and treat 
each subspace independently. 
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3.3 Convergence 


We briefly discuss convergence of the solution of the discrete problem to that 
of the continuum problem. We assume that (55) holds. Inserting the exact 
smooth solution u(t,x ) into the scheme v n+l = Qv n generates truncation 
errors as inhomogeneous terms in the difference approximation and in the 
initial data. The error grid-function ic” = u” — u(t n ,Xj) satisfies 


w] +1 = Qwj + F" , (56) 

= fj , (57) 

where F'‘ = (f>(t n , Xj)0(k Pl + h P2 ), and f J = ij)(xj)0(h P3 ) with <p smooth. 
The temporal accuracy of the scheme is p\ and the spatial accuracy is p 2 - The 
discrete version of Duhamel’s principle (see Theorem 5.1.1 in [20]) gives the 
estimate 


\\w n \\ h , D+ <Ke at " (j|™°|| M3+ + II^IImn) < 0(k p1 + h p2 ), (58) 

provided that the initial data satisfies ||m 0 ||/ l! £) + < 0(h p2 ). If ^ is smooth this 
condition is satisfied and, in particular, it is satisfied for exact initial data. 

Inequality (58) implies convergence with respect to the discrete L 2 norm, 
lli^lU < despite the scheme being unstable with respect to this 

norm. Note that p-th order convergence is obtained, with p = min(p 1 ,p 2 ) 
assuming k = A h, even though the norm contains first order accurate one¬ 
sided difference operators. 


4 Applications 


In the following subsections we apply the theoretical tools discussed in Sec. 3 to 
different systems. We start with a first order strongly hyperbolic system with 
no lower order terms. We then investigate three second order in space systems: 
the wave equation, a generalization of the KWB formulation of Maxwell’s 
equations and the NOR formulation of Einstein’s equations. We show that 
the clear correspondence between strong hyperbolicity and the existence of a 
discrete symmetrizer which occurs in first order systems with no lower order 
terms is lost when the standard discretization is used for second order in space 
systems. Similarly, the simple correspondence between characteristic speeds 
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and the von Neumann condition, Eq. (63), does not hold for second order in 
space systems. It is convenient to define the following quantities, 


2 • q £* 2 • 2 c f) ^X2 

x 0 = sm n * x =2^ sm &, n = 


q / -j o ’ 

i=i z 


i =1 


/l ’ 


(59) 


Note that the maximum of \q and X is Vd. We also recall that when the 
eigenvalues of P are imaginary, 


cr(kP) < a 0 a(Q) < 1, 


(60) 


where cio = 2 for ICN, a/8 for 4RK and \/3 for 3RK. 


4-1 Stability of first order strongly hyperbolic systems 


Our Erst application is a constant coefficient first order system in d spatial 
dimensions 


= y 'Ai&L 

dt dx i ’ 


(61) 


where u is a vector valued function of the space-time coordinates. We assume 
that the system is strongly hyperbolic and that it admits a symmetrizer, 
i.e., there exists a matrix H(u) in Fourier space, such that H(u)P(iu) + 
P*(iuo)H(tjj) = 0, where P{iu ) = i uJ l A l . The discrete symbol associated 
with the standard second order accurate discretization of this system is 

A(0 = T 5/ A% sin ti = PfilhP 1 sin£) , 

n i=1 


where we attached the subscript h to the discrete symbol to distinguish it 
from that of the continuum. We now construct the discrete symmetrizer 

= H{h~ l sin£). (62) 


Conditions (29)-(30) are satisfied and condition (28) is sufficient for stability. 
The latter becomes afikP) = Xxcr(A(n)) < a 0 , where A(n) = J2i=i n iA\ 
Hi = x' 1 sin , so that X/=i = 1- Since this inequality must hold for all /, 
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and the quantity x reaches its maximum value \fd at £* = ±vr/2, we obtain 
the stability condition 


A < 


ot o 

cr(A(n))y/d 


(63) 


In the symmetrizable hyperbolic case one can take the discrete symmetrizer 
to be the same as that of the continuum (which, by definition, is independent 
of u ) 

MO = H. (64) 


This analysis of first order strongly hyperbolic systems shows that if the char¬ 
acteristic speeds depend neither on the direction nor on the dimensionality of 
the problem, i.e., if <j(A(n)) depends neither on n nor on d, then the Courant 
limit has a 1 /\fd dependence. In addition, when the second order accurate 
centered difference operator D 0 is used to approximate the spatial derivatives, 
a Courant limit violation would manifest itself as a rapid growth of the mid 
high frequency mode |£j| = f ~ 1.571. This mode is present if N is a mul¬ 
tiple of 4. A similar analysis shows that in the fourth order accurate case 
the situation differs. The Courant limit is 1.372 times smaller than (63) and 
above this limit the most rapid growth occurs at a slightly higher frequency, 
|£i| = 2arctan(6 1//2 /(4 — 6 1 / 2 )) 1 / 2 ~ 1.797. See also Appendix B. 


4-2 First order in time and second order in space wave equation 


In this section we discuss the stability properties of an approximation of the 
d dimensional wave equation written as a first order in time and second order 
in space system 


dt<j>(t,x) = U(t,x ), 

(65) 

<9 f n (t,x) = J2di(f)(t,x). 

(66) 


i=1 


In the introduction we pointed out that the Cauchy problem for this system 
is not well-posed in L 2 . One can expect that a direct application of definition 
(12), which is based on the discrete L 2 norm, to a scheme approximating (65)- 
(66) would lead to the conclusion that the scheme is unstable. The first order 
reduction, however, is well-posed in L 2 (it is symmetric hyperbolic), hence the 
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second order system satisfies an energy estimate with respect to 
\4>(x,t)\ 2 + \U(x,t)\ 2 + J2 \di(p{x,t)\ 2 d d x . 

i=l 


wv)ir=/ 


( 67 ) 


In this section we show stability for the standard discretization of this system, 
both by the pseudo-discrete reduction method given in Sec. 3.2, and by a 
direct discrete reduction in physical space. The two methods give equivalent 
results. 

Following the method of lines, we first discretize space and leave time contin¬ 
uous, 


^fe(i) = n (68) 

4n ; ,'(i) = EB +i B-^(i). (69) 

ai j=l 


Using the technique described in Sec. 3.2, we see that the (principal) symbol 
of the second order semi-discrete problem 


1 0 1 N 

T -i = 

( in i\ 

K -n 2 o y 


^ —ifl 1 J 


(70) 


has purely imaginary eigenvalues ±zf2. The matrix T diagonalizes P. The sum 
of the squared moduli of the characteristic variables gives the conserved energy 
(here D = 1/2J) 


v* (T -1 )* DT~ x v = \i£l(j> I 2 + |ri| 2 = n 2 \4>\ 2 + |n| 2 . (71) 

By taking K = 1 in (53) we see that we have numerical stability with respect 
to the discrete norm 

IMHd + = + n| + ’t(D +i h) 2 )h d , (72) 

j i= 1 

provided that the von Neumann condition 

A < a 0 /(2Vd) , (73) 

which follows from a(kP) = kfl = 2A%2 < ckq, is satisfied. 
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We now illustrate a different method for proving stability of this system. A 
discrete reduction to first order can be performed before going to Fourier space. 
We introduce the quantities 

Xf = D+ifr (74) 


and obtain the reduced system 


J t m = n,(i), 

jXf{t) = D +< TL,(t). 


(75) 

(76) 

(77) 


Notice that only if Eq. (74) is identically satisfied is the reduced system 
equivalent to the original one. It is important to check whether the evolu¬ 
tion equations (75)-(77) are compatible with this requirement. Let C^\t) = 

Xj 1 ’ — If we prescribe initial data such that C^\ 0) = 0, then at later 

times Cj l \t ) = 0. This is a consequence of the fact that 

= 4(xf (t) - D+*(i)) = 0 . (78) 


There is a one-to-one correspondence between solutions of (68)~(69) and those 
of (74)-(77). Furthermore, one can check that the time integrator does not 
spoil the propagation of the constraints. 

Ignoring lower order terms, the symbol associated with the reduced system 
(75)-(77) is anti-Hermitian, therefore Eq. (30) is satisfied with H = 1. The 
non-trivial eigenvalues of P are ±zf2, the same as those of the original system 
(68)-(69). This proves that (73) is a necessary and sufficient condition for 
stability with respect to the discrete norm (72). 

This specific discrete reduction to first order, and the pseudo-discrete reduc¬ 
tion to first order described in Sec. 3.2 give equivalent results. 


4-2.1 Fourth order accuracy 

In hyperbolic problems a fourth order accurate spatial discretization requires 
significantly fewer grid-points per wavelength for a given tolerance error (see 
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[20] and appendix B). The stability proof for the fourth order accurate dis¬ 
cretization of the d-dimensional wave equation 


= (79) 

2 n ,.(i) = y D +i D_ t (1 - y 2 d *‘ d -<) (80) 


is similar to the second order accurate case. The discrete symbol and diago¬ 
nalizing matrix are 


P 


' o k 

T -l = 

(,A l\ 



\~ iA V 


(81) 


where A 2 = ^ Yli=i sin 2 §■ (l + § sin 2 §■) > has purely imaginary eigenvalues 
±iA. Taking D = 1/2/ we get the conserved quantity 

(T” 1 u)*DT _ 1 f) = A 2 |0 | 2 + |fl| 2 . (82) 


Since fl 2 < A 2 < |D 2 , by taking K = 4/3 in (53) we see that we have 
numerical stability with respect to the norm (72) provided that the principal 
symbol P satisfies a(kP) < a 0 - This gives a stability limit of A < y/3a 0 / (Ay/d). 


4-2.2 A note about the Do-norm and the Dq discretization 


Replacing the one sided difference operators D + i with centered difference op¬ 
erators Dot in the norm (72) leads to difficulties, as the Do-norm does not 
capture the highest frequency mode. In fact, it is possible to construct a family 
of solutions of (68)-(69) proportional to (—1) J for which the D 0 -energy esti¬ 
mate fails. For this purpose it is sufficient to consider <f>j(t) = (—l)- 7 cos(2 t/h), 
Ilj(t) = — 2/h(— l )- 7 sin(2£//i), which gives 


l^(*)IU,r > 0 

^(0) |U, D 0 



2 1 4 „ 

A + i? Sin ” 


21\ 1/2 

h 


(83) 


where ||v(0llft,D o = +A 2 - + (D 0 (pj) 2 )h. It it not possible to find constants 

K and a such that the ratio is bounded by Ke at , independently of the space 
step h. 

It has been suggested that the use of Dq rather than D + D_ for the second 
spatial derivatives may improve the stability properties of a second order in 
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space scheme [24,25]. To investigate this we study the wave equation in one 
space dimension discretized as 



(84) 

jll = 

(85) 


The eigenvalues of kP are ±?'A sin £, which shows that the von Neumann con¬ 
dition is satisfied as long as A < «o- Both the stencil and the maximum time 
step compatible with the von Neumann condition are twice what they are for 
the D + D. _ discretization. However, for a given spatial resolution the numerical 
speed of propagation has an error which is four times that of the D + D_ case 
(see Appendix B). 


So far, we have only shown that the scheme is unstable if A > ck 0 - By looking 
at the discrete symbol 


1 0 l\ 

V -A s in 2 £ 0 J 


( 86 ) 


we see that there might be a problem for |£| = tt. In this case the symbol is 
not diagonalizable. To explicitly show that the system (84)-(85) is unstable 
with respect to the norm 

IMI l,D+ = E (fi + n ? + P+E 2 ) h ( 87 ) 

j 


it is sufficient to consider the family of initial data (f>j(0) = 0,^(0) = (—l)- 7 ’, 
generating the solution (f)j{t) = = (—1) J . As h —► 0 the ratio 


ll^( t )IU,g+ = A , f 2 , 1/2 

||w(0)|| h) r» + V h2 J 


( 88 ) 


grows without bound. 

Had we chosen the £) 0 -norm, however, we would have concluded that the 
scheme satisfies the required estimate. This is because this norm does not 
capture the highest frequency mode (f)j = (—1) J . A desirable requirement of 
a norm is that if a scheme is stable with respect to that norm, then it will 
remain stable with respect to the same norm when perturbed with lower order 
terms (independently of how these are discretized). The modified problem 
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(89) 

(90) 


jn j (t) = D^ j (t) - D+<t>j(t) 

admits the family of exponentially growing solutions (f>j(t) = (-l)- 7 exp(^/2/ ht ), 
n ? (f) = (—iy^j2/hexp(^j2/ht) which leads to unbounded growth in the ratio 


M*)lkr> 0 

^(0) |U,Z?o 



(91) 


If we want to be able to decide whether a scheme is stable or not just by 
looking at the principal part of the discrete system, then we must conclude 
that the Do-energy is not a suitable energy. 

We note that the requirement that stability should not depend on how lower 
order terms are discretized was crucial. If we restrict ourselves to the pertur¬ 
bation D 0 (ftj, then the scheme is still stable with respect to the D 0 -energy. If 
one wants to be able to discretize lower order terms freely, as we do, then one 
is forced to reject the Dq discretization. 

Clearly it is the presence of high frequency modes that makes the Dq discretiza¬ 
tion unstable with respect to the D + -norm. The introduction of a mechanism 
that damps high frequency modes, such as artificial dissipation, may restore 
stability. In the system 


j t <t>j = - ah 3 (D+D_) 2 0j , 

jn, = D 2 „^-ah 3 (D + D_) 2 n, 


the same family of initial data used to prove instability of (84)—(85) gives 
||'y(i)||h,D + /||w(0)||h,r> + = (1 + t 2 + 4 1 2 /h 2 y^e~ 16at ^ h , which does not grow 
without bound. 


4-3 The generalized Knapp-Walker-Baumgarte system 


We now investigate more complex systems. We adopt the Einstein summation 
convention. We consider the KWB formulation of Maxwell’s equations [16] 


dtAi — —Ei , 


(92) 
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d t E t = -d k d k A t + d t T , 

<%r = o, 


(93) 

(94) 


and generalize it by introducing G = Y — rd k A k , giving 

d t Ai = -Ei , (95) 

d t Ei = -d k d k Ai + rdid k A k + d t G , (96) 

d t G = rd k E k . (97) 

For r = 0 we recover (92)-(94) and for r = 1 we obtain the Z1 system [17], 
which was recently introduced as a toy model for the Z4 formulation of General 
Relativity (see Sec. 4.6). We will show that although the parameter r plays 
no role at the continuum, at the discrete level it can have a severe impact on 
the stability properties. 


4-3.1 Continuum analysis 

If we Fourier transform (95)-(97) and introduce T = G + riuj k A k in place of 
G the system simplifies to 


dtAi — —Ei , 
d t Ei = uj 2 Ai + iuit , 

9tf = o. 

The eigenvalues and characteristic variables of the symbol are 


0, w (0) = f, 

±iuj, wl ± ' 1 — Ei T iajAi ± cDjf , 

where tD* = Ui/u and a r = J2 k =i u k- Note that the eigenvalues of the symbol 
are independent of the parameter r. To construct a conserved energy we take 
the combination 


Ec^llw^ + llw^ + alw^l 2 . 


To keep the notation compact we omit the sums. We need to check that this 
conserved quantity is equivalent to 8 

\u\ 2 =\E i \ 2 + oj 2 \A l \ 2 +\G\ 2 . 

8 From the results in Section 3 we only need to show that H is equivalent to I u , 
see inequality (36), which in this case means that there is no |^h| 2 term. 
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Since 


E c = | Ei\ 2 + (1 + a)|r| 2 + a; 2 |A| 2 - 2Re , 


we get 


\E{\ 2 + (1 + a — £i)|f | 2 + ^1 — — ^ uj‘ 2 \A i \'' < Ec 
— I^i| 2 + (i + a + £ A |r 1 2 + A 2 I Ai | 2 , 

where we used the inequality ±2Re(z!Z 2 ) < £|^i| 2 +£ _1 |^ 2 | 2 for £ > 0. Choosing 
a = 3/2, E\ = £2 1 — 2 gives 

K^\u\ 2 v < E c < Ki\ii\ 2 , 

with K\ = 3, where |u|p = \Ej\ 2 + oo 2 \Ai\- + |r| 2 . Using the inequality 

(1 - £)|^| 2 + (1 - £- 1 )|z 2 | 2 < |*1 + Z2? (98) 

< (1 + £)|^ l | 2 + (1 + £ 1 )|^ 2 | 2 , 

with £ > 0, we have that for any r, |w|p is equivalent to \u\ 2 , i.e. K 2 l |w|p < 
|ff| 2 < ih 2 |ix|f. We have the uniform estimate in Fourier space 

|*(f)| 2 < K 2 \u{t)\l < K 1 K 2 E c (t) = K X K 2 E C { 0) 

< K\K 2 \u (0 )\l < K 2 K 2 \u(Q)A (99) 

which implies the estimate in physical space with respect to the norm 

Hr-iiAir+ii^ir+RAir+iiGii 2 , (ioo) 

with no restrictions on the parameter r. 


4-3.2 Discrete analysis 

Consider now the semi-discrete system 

dt.Ai = —Ei , 

d t Ei = —D +k D_ k Ai + rD^A k + D 0i G , 
d t G = rD ok E k , 
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( 101 ) 

( 102 ) 

(103) 


(‘ 2 ') 

where D) k is the standard second order accurate approximation of the second 
partial derivative. The procedure is similar to that at the continuum. We 
Fourier transform and replace the variable G with Y = G + rj-sin^kAk and 
obtain 


dtAi — —Ei , 

d t Ei = ^ 0-(04 + , 

3*f = 0, 


where 0f (£) = Yll= l sin 2 % ~ r sin 4 

The eigenvalues of the matrix kP(£) and the corresponding characteristic vari¬ 
ables are 


0, w (0) = f , 

±2i0 i (e)A, 4 ±) = Ei T j®i{0Ai ± Si(0r, 

where 2Sj0,, = sin ^ ?; . The requirement that a(kP) < «o imposes the restriction 
r < 1 on the parameter. If this condition is violated, then the semi-discrete 
scheme is unstable (and the fully discrete scheme would be unconditionally 
unstable). Furthermore, for r = 1, which corresponds to the Z1 system, the 
matrix P(±7r, 0, 0) (corresponding to the highest frequency in the x direction) 
is not diagonalizable and one can show that the system admits frequency de¬ 
pendent linearly growing solutions which violate the discrete energy estimate. 

Assume r < 1. The expression 


E c = 


\w 


(+)|2 


+ 


\W 


(—) 12 


+ a\Y\ 


= \Ei\ 2 + (a + s 4 2 )|f| 2 + -^0-|A| 2 


2Re ( — sin ^AiY 


is conserved. We want to show that it is equivalent to \ii 

|Gf. 


E i \ 2 + n 2 \A i \ 2 + 


We first show that E c is equivalent to |w|p = \Ei\ 2 + Q-\Ai\ 2 + \Y\ 2 . We 
distinguish now between two possibilities: r < 0 and 0 < r < 1. In either case 
we have that |s*| < 1. In the first case, using the inequality X 2 — ©i — (1 — r )x\ 
we get 


— \Ei \ 2 + (a + 1 + £2)!!- h + — r + ~^ xl\A\ 2 - 

If we take, for example, a > 3, £1 = 2, e 2 = 1/2, then there exist constants 
K\ and K- 2 such that Ad | u | 2 < Ec < X 2 |w|r- 

For the case 0 < r < 1, using the inequality (1 — r)x\ < © 2 < xi we S e ^ 

\E{\ 2 + (a — £i)|I y + ^1 — r — — ^ xl\A \ 2 < < 

l-^/l 2 + (a +1 + s 2 ) |r | 2 + — ^1 + — ^ x 2 |A| 2 - 

If we choose a > e± > 1/(1 — r) we have the equivalence to |£t|p. On the other 
hand, using 

sinf fc | < |fi|, (104) 

one can show that the norms |u|p and Inti 2 are equivalent. This proves stability 
with respect to the norm 


(IIAIIS 


+ l|£ 


i\\h 


\D+ k A . 112 1 "^" 2 


'+ksn\\h 


+ II^IIj 


1/2 


(105) 


Note that the Cauchy problem for the continuum system is well-posed for all 
values of r, but the discrete system is stable only for r < 1. For r < 1/2 the 
von Neumann condition gives a Courant limit of A < «o/(2\/3 — r). Moreover, 
the numerical speeds of propagation depend on r. 


4-4 The Nagy-Ortiz-Reula system 


The NOR formulation of Einstein’s equations linearized about Minkowski 
space with zero shift and densitized lapse (a = det^j) 1 / 2 ) has the form 



(106) 

dtKij = ~d k d k Xij + '-didjT + dpfp , 

(107) 

d t fi = rdiK , 

(108) 


where r = S klr jki- This system corresponds to the one in [10] with the choice of 
parameters a — b — a — 1, c = 0 and p = r + 2. It is obtained from the ADM 
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system with densitized lapse by introducing the variables fi = dflij — diT, 
which are used in the evolution equations for the K t] variables, and adding 
the momentum constraint to the time derivative of the new variables. 


4-4-1 Continuum analysis 

We Fourier transform the system and introduce r, = /,; + obtaining 

dtjij — Crj , 

dtKij ^ Tp "F , 

d t fi = 0. 


The eigenvalues and characteristic variables associated with the symbol are 

0, w^ = ti, 

±icv , w\f } = Kij =F ± u)(jf j). 

Proceeding in the usual manner we construct a conserved quantity and show 
that it is equivalent to 

\ii\ 2 = \k i3 \ 2 + UJ 2 \%\ 2 + \f l \ 2 . 


We have 


E c= jK 
= I K 12 


(+)|2 


+ 2^ 


(—) 12 


+ a\w. 


( 0 ) ,2 


1 


ij i T |(jTy) | T |7ul ) T cz.|Ty i 


Since 


0 < |%fj)| 2 < \Citjl 2 < |fj| 2 - —|7p-| 2 - £i|f*| 2 < —2Re [iufiijtj 

£l \ 

<-|7of+ ^|f l | 2 , 

£'2 

we obtain the equivalence with |u|p, 


30 


I Kij | 2 + ^ — — ^ w 2 \%\ 2 + (a — £i)|fj| 2 < Eq 

— \k-ij\~ + ^ + (i + a + ^2) ir^ | 2 , 

by choosing a = 3, £1 = 2, e 2 = 1. Finally, noting that \t | 2 < 3|7y| 2 one can 
show that \u\y and |w| 2 are equivalent. 


4-4-2 Discrete analysis 

We consider the standard second order accurate discretization of system (106)- 
(108). The semi-discrete system is 

d tlij = - 2Kij, (109) 

dtKij = ~~D +k D_ k r j i j + -D\^t + Dofifj ), (HO) 

d t f% — rD 0i K . (Ill) 


Taking the Fourier transform and introducing = fi + sin £,;f gives 


dtiij 2 Kjj , 

T * % a. 

dtKij 9 ^ T ~\~ — sin£(jFj), 

<9/f* = 0, 

where 

a f 0 * ^ •? 

4 . 4 6 . . ' 

l-^ sm 2 ' = * 


The eigenvalues of kP and the corresponding characteristic variables are 


0, 

.(0) 

w\ 

±2i0A, 


±2fy 2 A, 





= Ti 


1 

t 


sin jj 

20 


r®, 


- 1 , sin^r,-) 

= Kij =F —ikl'-fij ± - - —— 

J 2 2 X 2 


, 3 , 

TF 
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where 0 2 = xl ~ r ELi a h = sin4 f, °t K u = K«, °ilu = 7»> <R 4r i = f 1 *, 
and = (A? - SijA/3). 

Note that stability demands that r < 1 (p < 3). Furthermore, the von Neu¬ 
mann condition depends on the value of this parameter. Explicitly, this is 


2 max|^|<^{0, X 2 } 


and its dependence on r is illustrated in Figure 1. This is in contrast to the 
fact that at the continuum r has no influence on the characteristic speeds or 
the hyperbolicity of the system. 



Fig. 1. The von Neumann condition for the second order accurate discretization of 
the NOR system in 3D using 4RK as a function of the parameter r. For r > 1 the 
scheme is unconditionally unstable. 

We now restrict ourselves to the case r = 0 and prove numerical stability. In 
this case the characteristic variables associated with the non trivial eigenvalues 
are 


w-f = Kij T 2 ^% ± 


sin^rp 

2y 2 


( 112 ) 


A conserved quantity is 


E c = -| w. 


(+)|2 

ij I 

2 


+ 2 lW ‘> 
2 


(—) 12 


( 0 ) |2 


+ a | tu¬ 
ft 2 2 (l 


— \Kij\ + I s (iFy)| + ~\lij\ ~ Re ( — sin^y^Tj ] + a|Fj| , 


where 2y 2 = sin^. 


Since 
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N < i 

0 < |s (i f j} | 2 < ISifjl 2 < |f ,;| 2 - ^2Xl\lij\ 2 - £l|r. t | 2 
< —2Re ^sin^f^ < + £ 2 |fi| 2 , 

we have the equivalence with |w|p. Inequality (104) guarantees the equivalence 
of the latter with This completes the proof of stability with respect to 
the norm 


+ \\K, 


1-3 


I \D+ k % 


13 


I 2 + 

\h + 


WfiWh) 


1/2 


(113) 


4-5 The ADM system 


With a densitized lapse function, a = det^^) 1 / 2 , the ADM equations lin¬ 
earized about the Minkowski solution in Cartesian coordinates take the form 


dtjij = -2Kij , (114) 

d t K i:j = d k d(ij j)k - ^d k d k -f - dfijT . (115) 

The symbol P{iuj ) of (114)-(115) is not diagonalizable and neither is that of 
its differential nor its pseudo-differential reduction. The family of solutions 
in which the only non vanishing components are 71 a = sin(u;x)i, K\a = 
— sin(tt;x)/ 2 , where A = 2,3, can be used to explicitly show instability. It 
gives 


Mv)ll 

«(o, Oil 


(l + 4f 2 + 4u> 2 f 2 )' /2 , 


(116) 


where \\u(t, -)|| 2 = ||y^(t, -)|| 2 + || K^t, -)|| + \\d k ^ij(t, -)|| 2 . The ratio cannot be 
bounded by Ke at with K and a independent of to. 

To see that the second order accurate standard discretization is unstable we 
take 71 a = (—1 )H and K\a = (—l) J+ 1 / 2 . As in the continuum, the ratio 


l^(*)IU,r>+ 

^(0)||h,D + 


1 + 4t 2 + 16 


-) 

h 2 ) 


1/2 


(117) 
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cannot be bounded. We can nevertheless compute the von Neumann condition, 
which is given by 


A < 


a/3q;o 

2y/7d' 


(118) 


In [26] stability tests were done with the non linear version of this formulation. 
The domain used consisted of a thin channel, with an even number N of grid 
points in one spatial direction and 3 grid points in the other two directions. By 
taking this into account we see that modes corresponding to the frequencies 
£i = 7 r, and £2 = £3 = 27t/ 3 grow exponentially if A > 0.4163. Figure 2 in [26] 
confirms that with a Courant factor of A = 0.5 there is a violation of the von 
Neumann condition. 9 

Although the symbol associated with the continuum system (114) and (115) 
has four Jordan blocks of size two for any c j, interestingly, the symbol as¬ 
sociated with the semi-discrete problem obtained with the standard second 
order accurate discretization can have rather different properties. For Fourier 
modes traveling in directions parallel to the axis the continuum result still 
holds. However, for Fourier modes not parallel to any of the axis, we found 
that the symbol may have fewer Jordan blocks. For some Fourier frequencies 
we even noticed that the symbol is diagonalizable. There is no conflict be¬ 
tween this observation and the fact that the continuum problem is ill-posed. 
As shown at the beginning of this subsection the discrete initial value problem 
for the ADM system is also ill-posed. In the limit of high resolution, h —> 0 
(£ —> 0 and c v fixed), the discrete symbol is a perturbation of the continuum 
one 10 

P d = P c + 0(h 2 ). 


Even though for some frequencies P d is diagonalizable, the characteristic vari¬ 
ables become degenerate in the limit h —> 0 , which implies that the discrete 
symmetrizer becomes unbounded (it is not possible to find a K , independent 
of h , satisfying inequality (53)). 


9 A one-dimensional von Neumann analysis gives the limit (118) with d = 1 and 
a 0 = 2) which corresponds to 0.655. However, this would not capture the fact that 
there could be exponentially growing modes with non trivial dependence in the two 
thin directions. 

10 Note that in general by perturbing a non diagonalizable matrix one obtains a 
diagonalizable matrix, so the diagonalizability of the discrete ADM symbol for some 
frequencies should not be so surprising. 
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4-6 The Z4 system 


The same family of solutions that was used to show instability of the dis¬ 
cretized ADM equations can be used for the standard discretization of the 
linearized Z4 system [19] 


d t a = -f(K - ra0), 

dtKij didjOi -dkdk^ij ~\~ dk^(% r )j)k ~ 0^0jt ~h 2 d^Zj^ , 

9,e = i (a t a nu -dAr) + a k z k , 

dtZi = d k K lk - d,I\ + [),(-), 

for any values of the parameters / and m. This instability, however, is not 
present if the Dg discretization is used as in [24], in conjunction with the 
Do-norm. Furthermore, it is possible that artificial dissipation may cure this 
instability of the standard discretization, at least for 0 </^lorl = / = 
m/2, since in this case the continuum Cauchy problem is well-posed. Note that 
while we use the same family of solutions that was used to show instability for 
the ADM case, the two cases are very different: While the ADM instability is 
due to the lack of well-posedness of the continuum equations, the problem with 
the Z4 system arises purely at the discrete level, and can be traced back to the 
difference in structure between the principal symbols of the pseudodifferential 
first order reductions of the continuum and discrete equations, see Eqs. (43) 
and (52). For second order in space systems diagonalizability of the discrete 
symbol is not implied by diagonalizability of the continuum symbol. 

The ADM and Z4 examples suggest a simple criterion that can be used to rule 
out certain schemes. Any first order in time, second order in space system of 
PDEs which gives rise to an ill-posed problem when the first order and mixed 
second order spatial derivatives are dropped will result in an unstable scheme 
if the standard discretization is used and no artificial dissipation is added. This 
is a consequence of the fact that grid modes with the highest frequency belong 
to the kernel of the D 0 operator. Although the Dg discretization gives stable 
schemes with respect to the D 0 -norm, provided that the continuum problem 
is well-posed, it suffers from the limitations described in section 4.2.2. 


5 Testing stability 

When dealing with variable coefficient or non linear problems it can be dif¬ 
ficult, if not impossible, to prove stability with respect to a certain norm. 
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Numerical experiments are often the only option. Given a discretization of 
the linear initial value problem ( 1 ) and ( 2 ), a stability test should be aimed 
at establishing the existence of the constants a and K , independent of the 
initial data and for all h < h 0 (and possibly k < XqK), by computing the ratio 
between a suitable discrete norm at time-step t n = nk and its initial value, 

II l) n || 

<Ke atn . (119) 

||u u || 

Although it is not possible to infer stability by examining a finite number 
of numerical experiments (one would have to explore the entire set h < h o 
that appears in the definition of stability), it is usually not difficult to spot a 
trend of behavior as the resolution is increased. To ensure that a wide range of 
frequencies is excited, random initial data can be used [27], as no smoothness 
assumptions are used in the definition of stability. 

In the examples of first order in time, second order in space hyperbolic systems 
for which we are able to determine stability, we use a norm which is the discrete 
version of the continuum one. The derivatives are approximated using the 
one-sided operators D + (or, equivalently, D_) rather than D 0 . For the NOR 
system, for example, we use the square root of the expression 

Il7ijllh+ \\ K ij\\h+ H \\ D +klij\\h + WfiWh- 

i,j =1 i,j=l k,i,j =1 i =1 


If, as we vary the initial data and the resolution, the experiments indicate 
that the constants a and K in (119) exist, then one would conclude that the 
scheme appears to be stable. If not, the scheme appears to be unstable. 

In the nonlinear case, if the problem has a sufficiently smooth solution uq, 
then to first approximation the error equation can be linearized about Uq and 
convergence follows if the linearized equation is stable (Sec. 5.5 in [20]). Es¬ 
tablishing stability experimentally using the linearized equations would not be 
very practical. However, convergence to a known exact solution can be tested 
directly and it avoids many complications. Rather than testing for stability, 
one could test convergence in a more demanding way: initial data can be cho¬ 
sen which is not smooth, but is accurate to the correct order in the appropriate 
norm. For instance, for the NOR system, one would use the square root of 

E ll<b«llS+ E \mM+ E IP+it^Olli + EllVillt 

i,j =1 i,j=l k,i,j=l i =1 

where 5v = v — Uo, and one could add random noise to the initial data with 
amplitude h p for the K i3 and /) variables and h p+1 for the 7 ^ variables. The 
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scheme is convergent around the solution uq if the D + -norm of the error at 
time T is of order h p . In particular, this implies that for a convergent scheme 
the discrete L 2 norm of the error is of order h p if the D + -norm of the initial 
error is of order h p . 

Finally, we note that the notion of robust stability introduced in [27] does not 
imply nor follows from the concept of numerical stability investigated in this 
paper. 


5.1 Numerical tests 


We have performed numerical tests to complement the analytical stability 
results of Sec. 4. 

For each run, the numerical grid has dimensions 50px 4 x 4, where p = 1, 2,4, 8 
parameterizes the resolution, and we impose periodic boimdary conditions. 
The coordinate domain is x,y,z E [—0.5, 0.5). The time integrator is RK4 
with Courant factor A = 0.5. We choose random noise of order unity as initial 
data (except for the Z4 tests, see below) so that many discrete Fourier modes 
are present in the initial data. Empirically, we find that using smooth initial 
data in the constant coefficient problems of this paper can make it difficult to 
observe an instability. This was also noticed in the nonlinear case in [28]. 

Figure 2 shows the results for the ADM system. The apparent trend is that 
as the resolution is increased, and higher frequency Fourier modes are present 
in the initial data, the ratio of the D + -norm of the solution to its initial value 
at any given time increases. It appears that there is no K, a such that this 
quantity can be bounded by a function Ke atn , and this indicates that the 
system is unstable. 



Fig. 2. Linearized ADM stability test 
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In Figure 3 we show the results of the stability test for the linearized NOR 
system. The results suggest that the ratio of the D + -norm of the solution to 
its initial value remains bounded, and hence that the system is stable. This 
reflects the analytic result that we proved in Sec. 4. 



Fig. 3. Linearized NOR stability test, r = 0 

Showing the instability of the Z4 system was more complicated. In this case, 
it was not sufficient to use random initial data of order unity in all variables. 
When this was attempted, the ratio of the _D + -norm to its initial value re¬ 
mained bounded. In order to numerically demonstrate the instability, we used 
knowledge of the exact solution that violates the estimate. Random data of 
order unity was given to the variables K 22 and W 33 and the remaining vari¬ 
ables were set to zero. The test results for the linearized Z4 system are shown 
in Figure 4, and confirm that this system is unstable. 



Fig. 4. Linearized Z4 stability test, / = l,m = 2. Initial data consists of random 
values in K 22 and -K 33 , and all other variables are zero. 

When artificial dissipation with a = 0.02 is used, the linearized Z4 system 
tested with the same initial data shows no sign of instability. See Figure 5. 
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Fig. 5. Linearized Z4 stability test with dissipation cr = 0.02, / = l,m = 2. Initial 
data consists of random values in K 22 and /L 33 , and all other variables are zero. 

The example of the Z4 system shows that numerical testing of stability is not 
always straightforward, and that schemes which appear stable for simple test 
cases may in fact be unstable. All tests were done using the standard second 
order accurate discretization. 


6 Discussion 


In this work we extended the notion of numerical stability of finite difference 
approximations to include hyperbolic systems that are first order in time and 
second order in space. We considered the standard discretization of the wave 
equation, a generalization of the KWB formulation of electromagnetism and 
the NOR formulation of Einstein’s equations linearized about the Minkowski 
solution. By analyzing the symbol of the second order system, and construct¬ 
ing a discrete symmetrizer, we were able to prove stability in a discrete norm 
containing one-sided difference operators, provided that the von Neumann 
condition is satisfied. Consistency and stability with respect to the D + -norm 
imply convergence with respect to the discrete L 2 norm. We also found that 
in some cases (r > 1 in the NOR and generalized KWB systems, and Z4) 
standard discretizations of well-posed continuum problems can lead to un¬ 
conditionally unstable schemes. This is closely related to the instability of the 
fully second order shifted wave equation investigated in [29], but our examples 
contain no shift terms. 

Our analysis of discretizations of first order in time hyperbolic systems shows 
that in the first order in space case there is a clear correspondence between 
strong hyperbolicity and numerical stability, and between characteristic speeds 
and Courant limits. See inequality (63) and Eq. (64). In the second order in 
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space case, on the other hand, the mixing of D± and D 0 operators breaks 
this correspondence. To restore the correspondence one could use the Dq dis¬ 
cretization, however, as discussed in Sec. 4.2.2, this can lead to difficulties. 

In Sec. 4.6 we propose a simple criterion that can be used to rule out certain 
schemes when the standard discretization is used and no artificial dissipation 
is added. This criterion detects schemes in which the highest frequency mode 
grows faster as the resolution is increased. 

We also discuss stability tests for second order in space systems. These tests 
should be aimed at establishing the existence, for sufficiently small h, of the 
constants K and a that appear in the definition of stability with respect to 
the Z} + -norm. In the nonlinear case the situation is more complicated. In this 
case we suggest, when an exact smooth solution of the continuum problem 
is available, to do convergence tests with initial data given by that of the 
continuum problem plus random noise of order h p with respect to the D + - 
norrn (see Sec. 5). 

Although our analysis was restricted to the constant coefficient case, we expect 
that for the variable coefficient case generalizations of results similar to those 
presented in Sec. 6.6 of [20] for first order hyperbolic systems, where artificial 
dissipation plays an important role, might apply. 
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A Time integrators 


In this work we restrict our attention to the following three time integrators: 
3rd and 4th order Runge-Kutta, and iterative Crank-Nicholson [30]. Given a 
system of ordinary differential equations, dy/dt = f(t,y(t)), these integrators 
are defined as 

3RK 


k l = kf{t ni y n ) 

k 2 = kf(t n + kj2 , y n + fci/2) 
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h = kf(t n + 3k/A, y n + 3/c 2 /4) 
y n+l — y n + (2/ci + 3k 2 + 4/c 3 )/9 

4RK 

k 1 = kf(t n ,y n ) 
k 2 = kf(t n + /c/2, + ki/2) 

k 3 = kf (t n + /c/2, + Zc 2 /2) 

h = kf(t n + k,y n + k 3 ) 

y n +l =y n + ( ki + 2h + 2h + h y 6 

ICN 

fci = kf(t n ,y n ) 

k 2 = kf(t n + /c/2, y n + ki/2) 

h — kf ( t n + /c/2, y n + k 2 / 2) 


B Some numerical properties of first and second order systems 


In this section we assume that the time integrator is one of those discussed 
in Appendix A. We consider standard second and fourth order accurate dis¬ 
cretizations of the following two toy model problems 

Ut u x j (B.l) 


and 


ft n, ir fxx ■ 


(B.2) 


Eq. (B.l) arises in the full reduction to first order of f tt = f xx , while (B.2) 
represents its reduction in time. If we denote by A(£) an eigenvalue of the 
discrete symbol, the corresponding phase and group velocities are given by 


v p = i 


,A(0 


UJ 

^ak), 


where £ = uh. In the following table we compute the numerical phase ve¬ 
locities, Up, group velocities, v g , the Courant limits (C.I.), the frequencies of 
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undamped modes (u.m.) and of the first unstable mode (f.u.m.) for the two 
systems. The numerical phase and group velocities are plotted in Figure B.l 
as a function of £. 

In the table we used A 2 = 1 + | sin 2 The exact continuum phase and group 
velocity is 1. The Taylor expansion of the numerical velocities gives an idea of 
the magnitude of the error, provided that enough grid-points per wave length 
are used. The table shows that in the second order accurate case the phase 
error for the wave equation is 4 times smaller than for the advective equation, 
and that this improvement in accuracy is even stronger for the fourth order 
accurate discretization. 


Furthermore, the standard discretizations of fully first order hyperbolic sys¬ 
tems have numerical phase velocities that vanish at the highest frequencies 
and numerical group velocities with the opposite sign to the continuum one. 
In numerical relativity simulations involving black holes which make use of 
the excision technique to handle the singularity one can expect to see nu¬ 
merical high frequency solutions escaping from the black hole, if a first order 
formulation combined with the standard discretization is used, unless artificial 
dissipation is added to the scheme. 


Finally, whereas for (B.l) the transition from second order accuracy to fourth 
order implies the reduction of the Courant limit by a factor of 1.372, for the 
second order in space system (B.2), this transition requires a Courant limit 
2 /a/ 3 ~ 1.155 times smaller. This indicates that there is an even higher gain 
in going to fourth order accuracy for second order in space formulations. 



2nd order accurate 

4th order accurate 


advective 

wave 

advective 

wave 

V P 

sin£ ^ 

£ ~ 

i-f + o(? 4 ) 

| sin | ss 

i - g + o(e>) 

±(i + ±>±) 

~ 1 - £ + o(e 6 ) 

| sin | A r; 

1 - ik + °(£ 6 ) 

v g 

cos£ ~ 

l-^ + OK 4 ) 

cos | ~ 

1-f + 0(? 4 ) 

1 — | sin 4 | « 

i-£ + o« s ) 

cos| (l ± | sin 2 /A 
« 1 - £ + °(£ 6 ) 

C.l. 

«0 

ao/2 

ao/1.372 

~ ao/2.309 

u.m. 

0,7r 

0 

0, 71 

0 

f.u.m. 

±f « ±1.571 

71 

±2 arctan ( f ,A ) 
W4-V6 / 

« ±1.797 

71 
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Fig. B.l. The phase (top) and group (bottom) velocities for the second (left) and 
fourth (right) order standard approximation of the advective equation (B.l) (red) 
and the wave equation (B.2) (green). 


C Discrete constraint propagation 


When simulating systems such as Maxwell’s or Einstein’s equations, one has 
to take into account that the data has to satisfy initial data constraints. The 
evolution equations guarantee that if these constraints are satisfied initially, 
then they will be satisfied at later times. In this appendix we show that even 
in the constant coefficient case, when using standard discretizations of second 
order in space systems, the discrete constraints do not propagate exactly. 
Initial data which satisfy the discrete constraints do not lead to constraint 
satisfying solutions. 
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As an example, we consider the ADM equations (114)-(115) with constraints 
C = ^ (didj^ij - didir) = 0 , Ci = djKij - d,K = 0 . 

For simplicity we confine ourselves to solutions which depend only on one 
space coordinate. The discretized constraints are 

C = —-D+D-^aa = 0 , C'i = —DqKaa = 0 , 

Ca = DqK]a = 0 , 

where A = 2,3. 

The time derivative of the first constraint cannot be expressed in terms of 
finite difference combinations of the constraints 

jC = d + d_k aa ± -AA. 

This is to be contrasted with the fact that in the constant coefficient case, 
the discrete constraints of a first order reduction would propagate as in the 
continuum, with partial derivatives replaced by D 0 operators. Furthermore, 
this issue would not be present if one used Dq to approximate the second 
derivatives. 
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